# The following analyses were carried out using R version 4.0.5 (2021-03-31)

# > sessionInfo()
#R version 4.0.5 (2021-03-31)
#Platform: x86_64-pc-linux-gnu (64-bit)
#Running under: Red Hat Enterprise Linux 8.9 (Ootpa)

#Matrix products: default
#BLAS/LAPACK: /opt/pub/eb/apps/cascadelake/OpenBLAS/0.3.12-GCC-10.2.0/lib/libopenblas_skylakexp-r0.3.12.so

#locale:
#[1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#[3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
#[5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
#[7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
#[9] LC_ADDRESS=C               LC_TELEPHONE=C            
#[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

#attached base packages:
#  [1] stats     graphics  grDevices utils     datasets  methods   base     

#loaded via a namespace (and not attached):
#  [1] compiler_4.0.5 tools_4.0.5   

###############################
## Install and load packages ##
###############################
library(dplyr)
library(ggplot2)
library(interplot)
library(grid)
library(gridExtra)
library(stargazer)
library(broom)
library(geepack)
library(lme4)
library(lmerTest)
library(Hmisc)
library(useful)
library(BBmisc)
library(magrittr)
library(MESS)
library(sjPlot)
library(ggeffects)
library(effects)
library(psych)
library(apaTables)


###############################
## Load Data ##
###############################
setwd("~/data")
reg.data1 <- read.csv('CalibratingAutonomy.csv')

###############################
## Calculating ICC ##
###############################
# Unconditional model: A model with no predictors
# Substantive question this model addresses: How much of the variation in perceived goverment quality is there within-and between- agencies?
katherine.unconditional.model <- lmer(Quality~1 + (1 | AGID),data = reg.data1) # This is an unconditional model that helps to evalaute if there is a significant proportion of varaince at agency level. If there is no significant varaince that can be attributed to cluster (i.e., agency), one simply can use a standard single-level multiple regression. In this model, we used the orginal data with missing values. In the subsequent analyses, we used compelete cases.  
summary(katherine.unconditional.model)

# ICC: Percentage of variance at agency level
Calculate_ICC <- function(my.lmer.object){
  variance <- my.lmer.object %>% VarCorr %>% as.data.frame
  ICC <- variance$vcov[1]/(variance$vcov[1]+variance$vcov[2])
  return(ICC) 
}                       # Function to calculate ICC. ICC can also be calculated by hand[variance at upper level divided by the sum of varaince at upper and lower levels]
Variance.at.agency.level<-Calculate_ICC(katherine.unconditional.model) %>%    # 
  round(2)
Variance.at.agency.level<- Variance.at.agency.level*100
Variance.at.agency.level


###############################
## GEE Models ##
###############################
# GEE model1: addressing the first research hypothesis without control variables 
reg.data1 <- reg.data1[order(reg.data1$AGID),] # We need to order by ID
reg.data1a<-reg.data1 # Rename data and standardize independence, capacity and discretion 
reg.data1a$Independence <- c(scale(reg.data1a$Independence))
reg.data1a$Discretion <- c(scale(reg.data1a$Discretion))
reg.data1a$Capacity<- c(scale(reg.data1a$Capacity))

reg.data1a$Education <- as.factor(reg.data1a$Education)

library(corrplot)
reg.data.corr  <- reg.data1a[,1:4]
apa.cor.table(reg.data.corr[,1:4], filename="CorrTableKeyVariables.doc") 

GEE.model1<-formula(Quality~Independence+ I(Independence^2) + Discretion +I(Discretion^2))
katherine.GEE.model1 <- geeglm(GEE.model1,id=AGID, corstr = "exchangeable", family = "gaussian", data=reg.data1a)# The correlation structure is assumed to be exchangable 
summary(katherine.GEE.model1)

# GEE model2: addressing the first research hypothesis by including control varaibles
GEE.model2<-formula(Quality~Independence + I(Independence^2) +Discretion +I(Discretion^2)+ Appointee+ Gender + Years + Education + Fed.District +Party) # # Independence^2" term was removed because was not significant
katherine.GEE.model2 <- geeglm(GEE.model2,id=AGID, corstr = "exchangeable", family = "gaussian", data=reg.data1a)
summary(katherine.GEE.model2)

# GEE resulsts tables: models 1 and 2 ( Main effects)

tab_model(katherine.GEE.model1, katherine.GEE.model2, dv.labels= c("GEE.Model1","GEE.Model2"), show.se = TRUE, file = "Table2.GEE.Models1en2.doc")

# H2:GEE
# Autonomy by capacity interaction 
GEE.model3<-formula(Quality~Independence*Capacity + Discretion*Capacity + Capacity*I(Discretion^2))
katherine.GEE.model3 <- geeglm(GEE.model3,id=AGID, corstr = "exchangeable", family = "gaussian", data=reg.data1a)
summary(katherine.GEE.model3)

# Autonomy by capacity interaction after taking the covariates into account 
GEE.model4<-formula(Quality~Independence*Capacity + Discretion*Capacity + Capacity*I(Discretion^2)+ Appointee+ Gender + Years + Education + Fed.District +Party ) #Independence^2" term was removed because was not significant in the previous model
katherine.GEE.model4 <- geeglm(GEE.model4,id=AGID, corstr = "exchangeable", family = "gaussian", data=reg.data1a)
summary(katherine.GEE.model4)

# Gee resulsts tables: models 3 and 4(Interaction effects)

tab_model(katherine.GEE.model3, katherine.GEE.model4, dv.labels= c("GEE.Model3","GEE.Model4"), show.se = TRUE, file = "Table2.GEE.InteractionModels3en4.doc")
# These two models mimic the lmer models tested above 

# Compare GEE models 
#Quasi-likelihood Information Criterion (QIC) is analogus to AIC. When comparing two GEE models, a model with lower QIC is preferred. As you can see here, perhaps removing some non-significant control variables is helpful (katherine.GEE.model1 is better than katherine.GEE.model2, same works for the interaction models(3 & 4)). 
mod1<-as.data.frame(QIC(katherine.GEE.model1))
mod2<-as.data.frame(QIC(katherine.GEE.model2))
mod3<-as.data.frame(QIC(katherine.GEE.model3)) 
mod4<-as.data.frame(QIC(katherine.GEE.model4))
QIC.Values<- cbind(mod1, mod2, mod3, mod4)
QIC.Values


# Interaction plots based on GEE model
p1a <- plot_model(model = katherine.GEE.model3, type = "pred", terms =  c( "Independence [all]" , "Capacity [meansd]"), axis.lim(y=c(1,5))) # "meansd" enables comparsion of  predictor varaibles at the mean,  1SD below, and 1SD bove od the moderator variable(.i.e., capacity)  

p1a <- p1a +
  theme_minimal() +
  theme(legend.position = "bottom",
        title = element_text(size = 8)) +
  ggtitle("Figure 3a: Impact of Indepedence on Quality \nat Varying Levels of Capacity") + 
  scale_colour_grey( # This change the fitted line and CI band to greyscale, you can adjust the color intensity by the start and end parameter
    start = 0.2,
    end = 0.8,
    na.value = "red",
    aesthetics = c("colour", "fill")
  )

###############################
## Figures ##
###############################

p1b <- plot_model(model = katherine.GEE.model3, type = "pred", terms =  c("Discretion[n=20]", "Capacity [meansd]"), axis.lim(y=c(1,5))) 

p1b <- p1b +
  theme_minimal() +
  theme(legend.position = "bottom",
        title = element_text(size = 8)) +
  ggtitle("Figure 3b: Impact of Discretion on Quality \n at Varying Levels of Capacity") +
  scale_colour_grey( # This change the fitted line and CI band to greyscale, you can adjust the color intensity by the start and end parameter
    start = 0.2,
    end = 0.8,
    na.value = "red",
    aesthetics = c("colour", "fill")
  )

grid.arrange(p1a,p1b, nrow = 1)
plots.sidebyside <- arrangeGrob(p1a,p1b, nrow = 1) 
ggsave(file="plot1.Interaction.GEE.png", plots.sidebyside)


# The complex plots: based on GEE
# The plot below is for GEE model 
p3 <- ggpredict(model = katherine.GEE.model3, terms = c("Discretion [n=20]", "Independence [-.5, 0, .5]", "Capacity [-1, 0, 1]"), hist=TRUE)
plot(p3) +
  theme_minimal() +
  ggtitle("") +
  scale_colour_grey( # This change the fitted line and CI band to greyscale, you can adjust the color intensity by the start and end parameter
    start = 0.2,
    end = 0.8,
    na.value = "red",
    aesthetics = c("colour", "fill")
  )
png(filename="Plot2.complex.GEE.png")
plot(p3) +
  theme_minimal() +
  ggtitle("") +
  scale_colour_grey( # This change the fitted line and CI band to greyscale, you can adjust the color intensity by the start and end parameter
    start = 0.2,
    end = 0.8,
    na.value = "red",
    aesthetics = c("colour", "fill")
  )
dev.off()



